##################################
# produce derivative sets from the original set 
# by picking a random fraction from the original
# and combining it with a random selection
# from the genome
##################################
derivativeSets<-function(inputGenes,numberSets,nGenesFromInput,nGenesFromGenome) {
    nGenesGenome=dim(inputGenes)[1];
    nGenesInput<-sum(inputGenes);
    
    # cannot take more than I have
    if(nGenesFromInput>nGenesInput) {
        nGenesFromInput<-nGenesInput;}
    
    #vector of all indices to take the random selection from
    #note: this could be improved by excluding the original genes from here
    allGenes<-c(1:nGenesGenome);
    
    # vector of input gene indices to take the sample from
    inputGeneIdxs<-which(inputGenes!=0);
    
    # initialize the matrix to hold the derivative sets
    derivSets<-rep(0,nGenesGenome*(numberSets+1));
    dim(derivSets)=c(nGenesGenome,(numberSets+1));
    
    # the first column is the original set
    derivSets[,1]<-inputGenes;
    # combine each derivative set
    for (i in 2:(numberSets+1)) {
            # pick from the original gene set
            sampleInput<-sample(inputGeneIdxs,nGenesFromInput);
            derivSets[sampleInput,i]<-1;
            sampleGenome<-sample(allGenes,nGenesFromGenome);
            derivSets[sampleGenome,i]<-1;
        }
    return(derivSets);
}

##################################
##################################
#--------------------------------------------------------------
# calculate overlap between two matrices 
# of binary gene/condition column vectors
# overlap is given as |intersect(a,b)|/sqrt(|a|*|b|)
# returns matrix of vector-vector overlaps 
#--------------------------------------------------------------
symmOverlap<-function(matrix1,matrix2) {    
    nGenesGenome<-dim(matrix1)[1];
    nVectors<-dim(matrix1)[2];
    # make binary
    matrix1[matrix1!=0]<-1;
    matrix2[matrix2!=0]<-1;

    # normalize each column in first matrix by 1/sqrt(nEntries)
    nEntries<-apply(matrix1,2,sum);
    # re-set zero size vectors to prevent div by zero   
    nEntries[nEntries==0]<-1;
    sizeMatrix<-rep(sqrt(nEntries),nGenesGenome);
    dim(sizeMatrix)=c(nVectors,nGenesGenome);
    sizeMatrix<-t(sizeMatrix);
    matrix1<-matrix1/sizeMatrix;

    # normalize each column in second matrix by 1/sqrt(nEntries)
    nEntries<-apply(matrix2,2,sum);
    # re-set zero size vectors to prevent div by zero
    nEntries[nEntries==0]<-1;
    sizeMatrix<-rep(sqrt(nEntries),nGenesGenome);
    dim(sizeMatrix)=c(nVectors,nGenesGenome);
    sizeMatrix<-t(sizeMatrix);
    matrix2<-matrix2/sizeMatrix;

    # now multiply normalized matrices to get overlap
    overlap<-t(matrix1)%*%matrix2;
    return(overlap);
}

##################################
# Threshold plot 
##################################
#--------------------------------------------------------------
threshplot<-function(similarityMatrix) {
    n1<-dim(similarityMatrix)[1];
    dim(similarityMatrix)=c(1,n1*n1);
    nPairings<-n1*(n1-1);
    thr<-rep(0,10);
    for (i in 1:100) {
        thr[i]<-(length(which(100*similarityMatrix>=(i-1)))-n1)/nPairings;
    }
    return(thr);
}

##################################
# Recurrent SigAlg 
##################################
#--------------------------------------------------------------
# 1) produce derivative sets
# 2) apply sigalg to each set
# 3) find matrix of overlaps between the output genes
recSigAlg<-function(inputSet,numberDerivSets,nGenesFromInput,nGenesFromGenome,threshC,threshG,nd1,nd2) {
    # 1) produce derivative sets
    derivSets<-derivativeSets(inputSet,numberDerivSets,nGenesFromInput,nGenesFromGenome);
    # 2) apply sigalg to each set
    output<-sigAlg(derivSets,nd1,nd2,threshC,threshG);
    # 3) find matrix of overlaps between the output genes
    outputGenes<-output[[1]];
    outputConds<-output[[2]];
    similarityMatrix<-symmOverlap(outputGenes,outputGenes);
    return(list(outputGenes,outputConds,similarityMatrix));
}
##################################
# Random control Recurrent SigAlg 
##################################
#--------------------------------------------------------------
recSigAlgRand<-function(inputSet,numberDerivSets,nGenesFromInput,nGenesFromGenome,threshC,threshG,nd1,nd2) {
    # 0) make a random input set
    nGenesInInput<-sum(inputSet);
    nGenesGenome<-dim(inputSet)[1];
    inputSet<-rep(0,nGenesGenome);
    dim(inputSet)=c(nGenesGenome,1);
    inputSet[sample(1:nGenesGenome,nGenesInInput)]<-1;
    # 1) produce derivative sets
    derivSets<-derivativeSets(inputSet,numberDerivSets,nGenesFromInput,nGenesFromGenome);
    # 2) apply sigalg to each set
    output<-sigAlg(derivSets,nd1,nd2,threshC,threshG);
    # 3) find matrix of overlaps between the output genes
    outputGenes<-output[[1]];
    similarityMatrix<-symmOverlap(outputGenes,outputGenes);
    return(similarityMatrix);
}


##################################
# Combined selected outputs 
##################################
#--------------------------------------------------------------
extractGenesAndConditions<-function(outGenes,outConditions,similarityMatrix,minRecc){
    # identify all output sets with recurrence >=minRecc
    checkVec<-100*similarityMatrix[1,];
    idxs<-which(checkVec>=minRecc);
    
    # consider only the recurrent output
    acceptedGenes<-outGenes[,idxs];
    acceptedConds<-outConditions[,idxs];
    
    # the output scores are averaged over the recurrent set scores
    geneScores<-apply(acceptedGenes,1,mean);
    condScores<-apply(acceptedConds,1,mean);
    return(list(geneScores,condScores));
}

##################################
# Apply Rec Sig Alg 
##################################
appRecSigAlg<-function(inputSet,numberDerivSets,fractionFromInput,fractionFromGenome,minRecc,threshC,threshG,nd1,nd2,fname) {	

    nGenesInInput<-sum(inputSet);
    nGenesFromInput<-round(fractionFromInput*nGenesInInput/100);
    nGenesFromGenome<-round(fractionFromGenome*nGenesInInput/100);
 
    # find derivative sets, apply signature alg. and find overlap matrices
    outReal<-recSigAlg(inputSet,numberDerivSets,nGenesFromInput,nGenesFromGenome,threshC,threshG,nd1,nd2) ;
    outRand<-recSigAlgRand(inputSet,numberDerivSets,nGenesFromInput,nGenesFromGenome,threshC,threshG,nd1,nd2) ;

    thrReal<-threshplot(outReal[[3]]);
    thrRand<-threshplot(outRand);

    bitmap(file=fname, type="png256", bg="transparent");
    plot(1:100,thrReal,type="b",col="blue",main="Recurrence plot",xlab="Similarity [%]",ylab="Fraction");
    points(1:100,thrRand,type="b",col="red");

    output<-extractGenesAndConditions(outReal[[1]],outReal[[2]],outReal[[3]],minRecc);
    return(output);
}